The Boson peak and the phonons in glasses 
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Abstract. Despite the presence of topological disorder, phonons seem to exist also in glasses at very high frequencies (THz) 
and they remarkably persist into the supercooled liquid. A universal feature of such a systems is the Boson peak, an excess 
of states over the standard Debye contribution at the vibrational density of states. Exploiting the euclidean random matrix 
theory of vibrations in amorphous systems, we show that this peak is the signature of a phase transition in the space of the 
stationary points of the energy, from a minima-dominated phase (with phonons) at low energy to a saddle-point dominated 
phase (without phonons). The theoretical predictions are checked by means of numeric simulations. 



INTRODUCTION 

X-ray and neutron scattering techniques allow to obtain 
very detailed physical insight into the high-frequency 
(0.1-lOTHz) vibrational dynamics of supercooled liq- 
uids and glasses. Within this range of frequencies their 
spectra reveal several universal properties |1|, related 
with the presence of sound-like excitations even for mo- 
menta p of the same order of magnitude of po, the first 
maximum of the static structure factor (typically corre- 
sponding to wave numbers of a few nm^^). This high- 
frequency sound is revealed as Brillouin-Uke peaks in the 
Thz region of the dynamic structure factor. An accessi- 
ble quantity to experiments is the vibrational density of 
states (VDOS), g{co), whose most striking feature is the 
presence of an excess of states over the Debye (O^ law 
in the "low" frequency region, (i.e. where the dispersion 
relation is linear, but still in the Thz region)|2|. This ex- 
cess of states is seen as a peak when plotting g{(o)/(0^ 
and has been named Boson peak (BP) ' . The peak posi- 
tion (Obp usually shifts to lower frequency on heating (2ll, 
except for the case of silica |4]. In this material the shift 
is seen on lowering the density |'5j . 

Due to its universality, the relevant physics underlying 
the Boson peak can be hopefully captured by some sim- 
ple model. Furthermore, several recent numerical simu- 



There exist alternative way.s of defining the boson peak from experi- 
ments, for example as a peak in Raman scattering data or as a peak in 
the difference between the observed VDOS of the glass and that of the 
corresponding crystal 



lations have shown that a model of harmonic vibrations 
is wholly adequate to describe this frequency range 1 6] 
and that anharmonicity need not be invoked. Given the 
presence of well formed local structures (Si02 tetrahe- 
dra, for instance) a natural approximation is to consider 
that the oscillation centers form a crystalline structure, 
the disorder in the atomic positions being mimicked by 
randomness in their interaction potential |7, ^ (disor- 
dered lattice models The main drawback of such a 
models is that they dramatically underestimate the scat- 
tering of sound waves 1 10]. A different approach studies 
vibrations around a topologically disordered |01 (liquid 
like) structure. It is followed by two different theories: 
modified mode-coupling theory 1 11 1 (which is not lim- 
ited to harmonic excitations), and euclidean random ma- 
trix theory (ERMT) 1 13 El El- ERMT owes its name 
to the fact that it formulates the vibrational problem as 
random matrix problem 1 15]. The matrices involved are 
called Euclidean random matrices ]16], and their study 
has required the development of new analytical tools. 
Both MCT and ERMT predict an enhanced scattering of 
sound waves as compared to disordered crystals. 

On the other hand, even within the harmonic frame- 
work the nature of the extra low frequency modes giving 
rise to the BP is still an open point. At a qualitative level, 
the frequency cObp is close to the loffe-Regel JItIi fre- 
quency (OiR, suggesting the possibility that the excess BP 
modes are localized ]18]. However, numerical simula- 
tions have shown that the localization edge is at frequen- 
cies greater than cobp and coir ] 19]. The loffe-Regel crite- 
rion signals rather a crossover to a region where the har- 
monic excitations are not longer propagating, due to the 



strong interaction with the disorder We call these modes 
glassons (since they do not propagate but "diffuse", they 
have also been called dijfusons 1 19|). A large bump of 
glassons is generally found around the loffe-Regel fre- 
quency, due to the flattening of the dispersion relation. 
This can be considered as the glass counterpart of the 
van Hove singularity of crystals |8, 20]. All the recently 
proposed theoretical frameworks predict that this peak 
of glassons should move to lower frequencies when ap- 
proaching an instability transition, where negative eigen- 
values (imaginary frequencies) appear. The aim of the 
paper is showing that the ERM theory makes sharp pre- 
dictions about the values of universal critical exponents 
describing the approach to this singularity and compar- 
ing them with numeric results. The emerging scenario 
describes the BP modes as given by the hybridization be- 
tween the phonons and the low-energy tail of the glasson 
peak which softens when the system approaches the in- 
stability transition i2lll22il . 



THE EUCLIDEAN RANDOM MATRIX 
THEORY 

The starting approximation is that particles can only os- 
cillate around fixed random positions, so that the position 
of particle / at time t is x, (f) = x"' + ^, (f ); the xj'' are 
quenched equilibrium positions (whose distribution must 
be specified) and ^,(f ) are the displacements. Hence the 
Hamiltonian is 
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H[x] = Y.Vix,-xj) c :,j:LM.^M^'>t^J (1) 

where the dynamical matrix M is an Euchdean Random 
Matrix: 



Mi, 



eg 



(2) 

with/fjv(x) = 5|,vV(x). 

In the one-excitation approximation the dynamic 
structure factor is 



Ep 



5(0)- toy,), (3) 



where e„ are the eigenvectors of the dynamical matrix 
and co„ its eigenfrequencies (= square root of eigenval- 
ues). The overline means average over the disordered 
quenched positions, whose distribution P[x'"'] has to be 
specified. The density of states (VDOS) is obtained in 
the limit of large momenta: 

g{co)=lim-P^S^'\p,co). (4) 



We can obtain 5'^'' (p, Co) through the resolvent G(p,z): 
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separating the axial tensor in a longitudinal term and a 
transversal one. The dynamic structure factor is then: 
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ImGL(p,a)2 + iO+). (6) 



A transverse dynamic structure factor (not measurable 
in experiments) can be defined in an analogous way. 
However, a most important and general result is that for 
p ^ °o the resolvent becomes isotropic: 
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= 5^v^Tr[z-M]-i. 



(7) 

So both longitudinal and transverse structure factors tend 
to a common limit (the VDOS, see eq. |4} at infinite 
momentum. ^ Leaving the potential y(r) unspecified and 
taking the simplified case P[x">] ^ 1 /V^ (V being the 
volume), one finds that: 
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_z-p/(0)+p/(p)-E(p,z) 
The self-energy E(p,z) is a matrix with the standard form 

E,v(p,z) =E,(;,)^+E,(;,) (V"^) • (9) 

which vanishes at p = oo and that can be computed in 
a series expansion in 1 /p . The main point is that the 
sum of all the infinite diagrams obtained recursively 
starting from this next-to-leading order result gives a 
self-consistent integral equation l22ll : 

^^v{p,z) ^^J ^^^V^A(q,P)GAa(q,^)^av(q,p). 

(10) 

where the vertices have the form V/j v (q, p) = P (//iv (q) — 
/|j V (p — q) ) ■ Let us remark that the self energy renormal- 
izes the dispersion relations and gives a finite width to the 
Brillouin peaks: 
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Hj) (/?)+ReEL,r(p,C0i,r(p)), 

\mT.Lj{p,(OL.T{p))l(OL,T{p)- (11) 



Conseguently, both the dispersion relations saturate at the same 
value. However due to the broadening of the line, they are rather ill- 
defined when CO ~ cOir 



The correlations between the equilibrium positions of 
the particles can be taken into account quite easily at the 
level of the superposition approximation in the above ap- 
proach. The results derived above for the case without 
correlations are translated to the correlated case by re- 
placing the functions /(x) by g^^'(x)/(x). In this way 
the usual power law divergence of the pair potential for 
|x| ^ is balanced by the exponential behaviour of the 
pair distribution function, and this ensures the existence 
of the Fourier transform of the product /(x)g'^) (x). 



The phase transition 

From equation (llOt it is possible to derive a few ana- 
lytic model-independent results about the arising of the 
Boson Peak. These results are expressed in form of scal- 
ing laws, whose exponents are predicted in this approx- 
imation. Simulations (see below) and experiments will 
allow to clarify the dependence of the exponents on the 
approximation. The VDOS can be obtained from 

g{(o) = -—lmG°°{(0^ + iQ+), (12) 
% 

where z = o)^ + /0+ and G°°(z) = Ump^ooG(p,z). 

Hence one have to solve the integral equation dlOl i in 
the p ^ °o limit: 

^ = z - p/(0) - pAG-{z) -pj P (q)G(q, z) 

(13) 

where G°°{z), A = {In)^^ f d^qf^{(i) and the last term 
are matrices proportional to the identity. 

The solution of the above integral equation yields a 
VDOS which contains both the phonons, since g{co) °^ 
0)^ at o 0, and the extended but not propagating glas- 
sons, described by a semicircle with center at © = p/(0) 
and radius 2y/pA 1 20]. If we Hmit (for pedagogical pur- 
poses) to the case where the VDOS changes because of 
changes in the density, the key point is the existence 
of a phase transition in the space of the eigenvalues of 
the Hessian matrix. In fact G°°(0) develops an imaginary 
part when p < p^ (Pc being a critical density), then the 
transition separates the stable phase (all positive eigen- 
values) and the unstable phase (negative and positive 
eigenvalues). The glassons are the modes which move 
towards the negative zone of the spectrum (hybridizing 
the phonons) when approaching the transition. The or- 
der parameter is (jO = — ImG°°(/0+) which vanishes as 
(p ~ lA]'^, with J3 = 1/2 and A = p - Pf 

The relation with the Boson Peak becomes quite clear 
when one writes down the VDOS in the stable phase 
arising from the theory without any reference to the 
control parameter, which then does not need to be p. In 



fact, one has 

,(«,A) = a.r„(a.A-P), H.)^[f2^ ^Jj , 

(14) 

with A defined in terms of an arbitrary control param- 
eter. The ERMT (in the cactus approximation) predicts 
p = 1, 7 = 3/2. Hence it exist a crossover frequency (in 
the region where the dispertion relation is still linear) be- 
tween a CO^ and a co^ region. We shall identify that with 
the BP frequency cobp- This implies that cobp ~ A^ and 
g{a)BP,A)/o)lp A"l, with 

i?=p(2-7)- (15) 

Let us note that the BP is indicated from a peak in the 
function g{(o)/(0^m not in g{(o). Summarizing, accord- 
ing to ERMT the BP frequency moves linearly toward 
when approaching the transition (from the stable side) 
and its height diverges as a power law whose exponent 
is T] = 1 /2. Eq. shows that at the level of the one- 
phonon approximation it can also be detected in the large 
p limit of the dynamic structure factor S{p,(o). 



BOSON PEAK IN A GAUSSIAN MODEL 

In order to confirm that the saddle-phonon transition de- 
scribed by the Euclidean Random Matrix theory is not 
an artifact of the approximation involved (cactus resum- 
mation), we solved numerically the cactus equation for 
the case where f{p) has a Gaussian form and compare 
with direct numerical results for the same model ll22ll . 
The model is described by 

Up) = Mp)^+Mp){s,.-^), 

krip) = i^-^j expi-p^/2alr). (16) 

This choice for /(p) is mainly due to its simplicity. How- 
ever the behaviour of the Boson peak close to the saddle- 
phonon transition have to be independent of the details 
of the model. Moreover the superposition approximation 
takes f^vip) = ■^[8{r)v^v{r)] yielding a fLjip) finite 
at p = 0, like in the Gaussian model. We shall consider 
various values of the density, which is here the control 
parameter, comparing the analytical (cactus) results with 
the numerical spectra and dynamic structure factor ob- 
tained from the method of moments i23ll . In the high 
density regime the approximations used in deriving the 
integral equation il3\ are quite good since the analytic 
solution reproduces the numerical spectrum (and in par- 
ticular the Debye behaviour) rather accurately (figGJ. 
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FIGURE 1. The VDOS g{0}) as a function of eigenfrequencies divided by the Debye behaviour co^ for p = 4 > p^, both numerical 
(obtained via the method of moments) and analytical. In the inset, we show g(a)) vs. CO. 



However, the crucial check regards the exponents of 
the transition. Figs.|2t and|2l; show that the position of 
the BP is linear with respect to A = (p — p^) and that 
the height of BP diverges as A^'/^. This confirms the 
theoretical predictions v = 1 and Tj = 1/2. In Fig. |2j5 
we determine the value of J by studying the fraction 
of unstable modes. In fact, in the region of parameters 
where p < p^ the fraction of unstable modes, defined as 
/» = ,/°oo^A (A)t/A, is given by 



/„(A)= rdcoco^gico/\A\)^ 
Jo 



|A| 



i+r 



(17) 



We find numerically (Figs. ^) that /„ ~ (p^ — p)^/^, 
i.e. 7=3/2. Finally, the order parameter (p vanishes as 
{Pc~p f with i3 = 1/2 (Fig. Ell). 

Hence our analytic treatment based on euclidean ran- 
dom matrix theory describe quite well the vibrational 
features of simple topologically disordered svstems ll22il . 
The following step is understanding to what degree of 
accurateness ERMT could describe the high frequency 
properties of more realistic systems [21J- 



BOSON PEAK IN A FRAGILE GLASS 

Starting from the hypothesis that the Thz region of su- 
percooled liquids and glasses can be described in terms 
of purely harmonic excitations, the origin of the Boson 
peak in glasses can be understood if we consider the 
ensemble of generalized inherent structures (GIS). For 
each equilibrium configuration the associated GIS is the 
nearest stationary point of the Hamiltonian. If we start 
from an equilibrium configuration at low temperature. 



the GIS is a local minimum, and it coincides with the 
more frequently used inherent structures (IS) |24](i.e. 
the nearest minimum of the Hamiltonian). On the con- 
trary, if we start from high temperature, the GISs are sad- 
dle points. In the GIS ensemble there is a sharp phase 
transition separating these two regimes. It takes place in 
glass-forming liquids |25| at the Mode Coupling tem- 
perature 1 26] (Tmc), above which liquid diffusion is no 
longer ruled by rare "activated" jumps between ISs but 
by the motion along the unstable directions of saddles. 
Phonons are present in the spectrum of the VDOS in 
the low temperature phase (IS dominated) but are ab- 
sent in the saddle phase. The key point is that the min- 
ima obtained starting from configurations below Tmc and 
the saddles obtained starting above Tmc join smoothly 
at Tmc- Thus we can study GIS as a single ensemble 
parametrized by their energy i25ll . 

Since this transition separates a phase where all the 
eigenvalues are positive from another one where even 
negative eigenvalues exist, we expect that ERMT is able 
to describe correctly this phenomenum. Hence we mea- 
sured numerically the values of some exponents pre- 
dicted by the theory in a simple model of a fragile 
glass ^2 LI. We simulated a soft-spheres binary mix- 
ture 12711 in the stable ^ghonon) phase with the Swap 
Monte Carlo algorithm |28], and computed the VDOS of 
the ISs obtained starting from equilibrium configurations 
at temperatures below Tmc^- 

In Fig|3lwe show that the theoretical predictions agree 
with the numerical data. Taking the IS's energy as the 



^ At very low T, where equilibrium is not achieved, runs were followed 
until e/5 got very close to its asymptotic value 
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FIGURE 2. Numeric results. The critical density in the following fits has been fixed to Pc = 0.54 and capital letters are the 
fitting parameters, (a): The position of Boson peak as a function of the density near the critical point, cogp vanishes linearly in 
A = p — Pf. (b): The fraction of unstable modes vanishes as (pc — p)^'^, thus yielding 7= 3/2. (c): The height of the BP, defined 
by g{cOBp) / CO^p, diverges as A^'', with T] ~ 1/2. (d): The order parameter (p = — ImG°°(0) vanishes as (pc — p)^, with jS --^ 1/2. 



relevant parameter for describing the spectral properties, 
one has A = Cc — e/5, e/5 being the energy of the ISs 
and gf the critical value. In fact, plotting g{(o)/(0^ a 
peak is clearly identified, which is seen to grow in height 
and shift to lower frequency on rising the IS's energy. 
Using all the spectra for which the peak position can be 
clearly identified, we find that the relationship between 
(Obp and the energy of the IS is linear (Fig. 3a). The 
energy at which (Obp becomes zero, e^, is found from 
a linear fit as = 1.74£ (£ is the energy scale), quite 
close to the value where the GIS stop to be minima (IS) 
and become saddles 1 23 1 . As for the height of the peak 
(Fig. 3b), the results are compatible with a power-law 
divergence. Fixing at the value 1.74e arising from 
the linear fit of cObp vs. ejs, a power-law fit yields an 
exponent /3 = 0.40(15), while if one fixes the exponent at 
/3 = 1 /2, then the critical value is ec = 1.752(2) e. Thus 
the numerical data are compatible with the theoretically 
predicted scaling, although we have not been able to 
work close to the critical point, and thus cannot get a 
great accuracy on the critical exponents or the critical 
point. 



CONCLUSIONS 

In summary, we have shown that the saddle (negative 
eigenvalues)-phonon (no negative eigenvalues) transi- 
tion, a common feauture of vibrating topologically disor- 
dered systems, is well described by the euclidean random 
matrix theory. It provides a coherent scenario for the aris- 
ing of a Boson Peak, which results from the hybridization 
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FIGURE 3. Scaling of the position, asp. and height of the 
Boson peak near the saddle-phonon transition (energies and 
frequencies in units of e and Wq respectively). Top (Obp is 
linear in the control parameter e/j and vanishes at ejs = Cc = 
1.74(1) e. Bottom The height of the Boson peak diverges as a 
power law with exponent /3 --^ 0.4. Height and position of the 
BP were obtained by fitting a parabola to the peak of g{(o) /aP'. 
Reprinted with permission from Nature k2h] . Copyright (2003) 
Macmillan Magazines Limited 

of acoustic modes with high-energy modes that soften 
upon approaching the transition. Hence we applied the 
theory to describe the saddle-phonon transition and the 
BP in supecooled liquids, comparing the predicted scal- 
ing laws with the numeric results obtained for a simple 
fragile glass former. The agreement found is quite en- 
couraging. The present discussion applies to experiments 
as long as one is in the regime where the inverse fre- 
quency is much larger than the structural relaxation time, 
when the harmonic approximation makes sense. We ex- 



pect that the saddle-phonon transition point of view will 
be able to bridge the realms of experiment and numeri- 
cal studies of the energy landscape. As a matter of fact, 
a recent experiment on the poly(methyl methacrylate) 
(PMMA) glass gave the first experimental confirmation 
of the ERMT predictions 
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